#Produce descriptive boxplots
boxplot(combined$PpoliticsNN~combined$CountryGr,data=combined, 
        xlab="Country group", ylab="Share of politics terms")
boxplot(combined$PpoliticsNN~combined$Crisis,data=combined, 
        xlab="Time relative to crisis", ylab="Share of politics terms")

boxplot(combined$PpoliticsNN~combined$Country,data=combined, 
        ylab="Share of politics terms", las=3)
boxplot(combined$PpolicyNN~combined$Country,data=combined, 
        ylab="Share of policy terms", las=3)
boxplot(combined$PeconNN~combined$Country,data=combined, 
        ylab="Share of economics terms", las=3)
boxplot(combined$PfiscalNN~combined$Country,data=combined, 
        ylab="Share of fiscal terms", las=3)
boxplot(combined$Pneutral~combined$Country,data=combined, 
        ylab="Share of neutral terms", las=3)


#Linear mixed effects models
install.packages("lme4")
library(lme4)
#Must run to get p values
library(lmerTest)
library(lmtest)

options(scipen = 999)
politics.full = lmer(PpoliticsNN ~ CountryGrR * Crisis +count+ (1|Country), data=combined)
summary(politics.full)
policy.full = lmer(PpolicyNN ~ CountryGrR * Crisis +count+ (1|Country) , data=combined)
summary(policy.full)
econ.full = lmer(PeconNN ~ CountryGrR * Crisis +count+ (1|Country) , data=combined)
summary(econ.full)
fiscal.full = lmer(PfiscalNN ~ CountryGrR * Crisis +count+ (1|Country) , data=combined)
summary(fiscal.full)

#Without random slopes, with unstandardized coefficients, without interaction
options(scipen = 999)
politics.noI = lmer(PpoliticsNN ~ CountryGrR + Crisis +count+ (1|Country), data=combined)
summary(politics.noI)
policy.noI = lmer(PpolicyNN ~ CountryGrR + Crisis +count+ (1|Country) , data=combined)
summary(policy.noI)
econ.noI = lmer(PeconNN ~ CountryGrR + Crisis +count+ (1|Country) , data=combined)
summary(econ.noI)
fiscal.noI = lmer(PfiscalNN ~ CountryGrR + Crisis +count+ (1|Country) , data=combined)
summary(fiscal.noI)

#Calculating R2
library(MuMIn)
r.squaredGLMM(politics.full)
r.squaredGLMM(policy.full)
r.squaredGLMM(econ.full)
r.squaredGLMM(fiscal.full)

r.squaredGLMM(politics.noI)
r.squaredGLMM(policy.noI)
r.squaredGLMM(econ.noI)
r.squaredGLMM(fiscal.noI)

#Plot effects
library(emmeans)
library(pbkrtest)

ref_grid(politics.full)
emmeans(politics.full, ~ CountryGrR | Crisis)
polplot<-emmip(politics.full, CountryGrR ~ Crisis, CIs = TRUE)
polplot<-polplot +scale_color_grey(start=0, end=0.6, name=" ")+theme_bw()+ylab("Expected share of terms related to politics")+xlab("Time period relative to crisis")
polplot

ref_grid(policy.full) 
emmeans(policy.full, ~ CountryGrR | Crisis)
policyplot<-emmip(policy.full, CountryGrR ~ Crisis, CIs = TRUE)
policyplot<-policyplot +scale_color_grey(start=0, end=0.6, name=" ")+theme_bw()+ylab("Expected share of terms related to policy")+xlab("Time period relative to crisis")
policyplot

ref_grid(econ.full) 
emmeans(econ.full, ~ CountryGrR | Crisis)
econplot<-emmip(econ.full, CountryGrR ~ Crisis, CIs = TRUE)
econplot<-econplot +scale_color_grey(start=0, end=0.6, name=" ")+theme_bw()+ylab("Expected share of terms related to economics")+xlab("Time period relative to crisis")
econplot

ref_grid(fiscal.full) 
emmeans(fiscal.full, ~ CountryGrR | Crisis)
fiscalplot<-emmip(fiscal.full, CountryGrR ~ Crisis, CIs = TRUE)
fiscalplot<-fiscalplot +scale_color_grey(start=0, end=0.6, name=" ")+theme_bw()+ylab("Expected share of terms related to fiscal category")+xlab("Time period relative to crisis")
fiscalplot


#Model not excluding neutral words
politics.sfull_wN = lmer(sPpolitics ~ CountryGrR * Crisis + (1|Country) , data=combined)
summary(politics.sfull_wN)
policy.sfull_wN = lmer(sPpolicy ~ CountryGrR * Crisis + (1|Country) , data=combined)
summary(policy.sfull_wN)
econ.sfull_wN = lmer(sPecon ~ CountryGrR * Crisis + (1|Country) , data=combined)
summary(econ.sfull_wN)
fiscal.sfull_wN = lmer(sPfiscal ~ CountryGrR * Crisis + (1|Country), data=combined)
summary(fiscal.sfull_wN)

#Model using absolute word count
politics.absfull = lmer(Spolitics ~ CountryGrR * Crisis + count + (1|Country) , data=combined)
summary(politics.absfull)
policy.absfull = lmer(Spolicy ~ CountryGrR * Crisis + count + (1|Country) , data=combined)
summary(policy.absfull)
econ.absfull = lmer(Secon ~ CountryGrR * Crisis + count + (1|Country) , data=combined)
summary(econ.absfull)
fiscal.absfull = lmer(Sfiscal ~ CountryGrR * Crisis + count + (1|Country), data=combined)
summary(fiscal.absfull)

#Eurozone analysis
combinedEUROZONE<-combined[ which(combined$Eurozone==1), ]
combinedEUROZONE$Core <- "Periphery"
combinedEUROZONE$Core[which(combinedEUROZONE$Country == "Austria")]<-"Core"
combinedEUROZONE$Core[which(combinedEUROZONE$Country == "Belgium")]<-"Core"
combinedEUROZONE$Core[which(combinedEUROZONE$Country == "Finland")]<-"Core"
combinedEUROZONE$Core[which(combinedEUROZONE$Country == "France")]<-"Core"
combinedEUROZONE$Core[which(combinedEUROZONE$Country == "Germany")]<-"Core"
combinedEUROZONE$Core[which(combinedEUROZONE$Country == "Luxemburg")]<-"Core"
combinedEUROZONE$Core[which(combinedEUROZONE$Country == "Netherlands")]<-"Core"
combinedEUROZONE$Core<-as.factor(combinedEUROZONE$Core)
summary(combinedEUROZONE$Core)

EUROpolitics.full = lmer(PpoliticsNN ~ Core * Crisis +count+ (1|Country), data=combinedEUROZONE)
summary(EUROpolitics.full)

ref_grid(EUROpolitics.sfull) 
emmeans(EUROpolitics.full, ~ Core | Crisis)
EUROpolplot<-emmip(EUROpolitics.full, Core ~ Crisis, CIs = TRUE)
EUROpolplot<-EUROpolplot +scale_color_grey(start=0, end=0.6, name=" ")+theme_bw()+ylab("Expected share of terms related to politics")+xlab("Time period relative to crisis")
EUROpolplot

EUROpolicy.full = lmer(PpolicyNN ~ Core * Crisis +count+ (1|Country) , data=combinedEUROZONE)
summary(EUROpolicy.full)
ref_grid(EUROpolicy.full)
EUROpolicyplot<-emmip(EUROpolicy.full, Core ~ Crisis, CIs = TRUE)
EUROpolicyplot<-EUROpolicyplot +scale_color_grey(start=0, end=0.6, name=" ")+theme_bw()+ylab("Expected share of terms related to policy")+xlab("Time period relative to crisis")
EUROpolicyplot


#correlations between DVs
cor(combined$PpoliticsNN,combined$PpolicyNN)
cor(combined$PpoliticsNN,combined$PeconNN)
cor(combined$PpolicyNN,combined$PeconNN)
cor(combined$PfiscalNN,combined$PeconNN)
cor(combined$PpoliticsNN,combined$PfiscalNN)
cor(combined$PpolicyNN,combined$PfiscalNN)

#Model diagnostics
plot(politics.full)
plot(policy.full)

qqnorm(residuals(politics.full))
qqnorm(residuals(policy.full))

ggplot(data.frame(lev=hatvalues(politics.full),pearson=residuals(politics.full,type="pearson")),
       aes(x=lev,y=pearson)) +
  geom_point() +
  theme_bw()
levId <- which(hatvalues(politics.full) >= .2)
combined[levId,c("PpoliticsNN","Country","Year","Report")]

ggplot(data.frame(lev=hatvalues(policy.full),pearson=residuals(policy.full,type="pearson")),
       aes(x=lev,y=pearson)) +
  geom_point() +
  theme_bw()
levId_policy <- which(hatvalues(policy.full) >= .2)


#Reveals Albania as highly influential for politics, but does not change substantitve findings
combinedNOAlb<-combined[ which(combined$Country!="Albania"), ]
politics.fullNOAlb = lmer(PpoliticsNN ~ CountryGrR * Crisis +count+ (1|Country), data=combinedNOAlb)
summary(politics.fullNOAlb)

ref_grid(politics.fullNOAlb) 
polplotnoRSnoA<-emmip(politics.fullNOAlb, CountryGrR ~ Crisis, CIs = TRUE)
polplotnoRSnoA<-polplotnoRSnoA +scale_color_grey(start=0, end=0.6, name=" ")+theme_bw()+ylab("Expected share of terms relative to mean")+xlab("Time period relative to crisis")
polplotnoRSnoA

plot(politics.fullNOAlb)
qqnorm(residuals(politics.fullNOAlb))

policy.fullNOAlb = lmer(PpolicyNN ~ CountryGrR * Crisis +scount+ (1|Country), data=combinedNOAlb)
summary(policy.fullNOAlb)
ref_grid(policy.fullNOAlb) 
policyplotnoRSnoA<-emmip(policy.sfullNOAlb, CountryGrR ~ Crisis, CIs = TRUE)
policyplotnoRSnoA<-policyplotnoRSnoA +scale_color_grey(start=0, end=0.6, name=" ")+theme_bw()+ylab("Expected share of terms relative to mean")+xlab("Time period relative to crisis")
policyplotnoRSnoA




